/******************************************************************************
 *
 * MantaFlow fluid solver framework
 * Copyright 2011 Tobias Pfaff, Nils Thuerey 
 *
 * This program is free software, distributed under the terms of the
 * Apache License, Version 2.0 
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Mesh edge collapse and subdivision
 *
 ******************************************************************************/

/******************************************************************************/
// Copyright note:
//
// These functions (C) Chris Wojtan
// Long-term goal is to unify with his split&merge codebase
//
/******************************************************************************/

#include "edgecollapse.h"
#include <queue>

using namespace std;

namespace Manta {

// 8-point butterfly subdivision scheme (as described by Brochu&Bridson 2009)
Vec3 ButterflySubdivision(Mesh& m, const Corner &ca, const Corner &cb)
{    
	Vec3 p = m.nodes(m.corners(ca.prev).node).pos + m.nodes(m.corners(ca.next).node).pos;
	Vec3 q = m.nodes(ca.node).pos + m.nodes(cb.node).pos;
	Vec3 r = m.nodes(m.corners(m.corners(ca.next).opposite).node).pos 
		   + m.nodes(m.corners(m.corners(ca.prev).opposite).node).pos 
		   + m.nodes(m.corners(m.corners(cb.next).opposite).node).pos 
		   + m.nodes(m.corners(m.corners(cb.prev).opposite).node).pos;
	return ( 8*p + 2*q - r)/16.0;
}

// Modified Butterfly Subdivision Scheme from:
// Interpolating Subdivision for Meshes with Arbitrary Topology
// Denis Zorin, Peter Schroder, and Wim Sweldens
// input the Corner that satisfies the following: 
//      c.prev.node is the extraordinary vertex, 
//      and c.next.node is the other vertex involved in the subdivision
Vec3 OneSidedButterflySubdivision(Mesh& m, const int valence, const Corner &c) {
	Vec3 out;
	Vec3 p0 = m.nodes(m.corners(c.prev).node).pos;
	Vec3 p1 = m.nodes(m.corners(c.next).node).pos;
	
	if(valence==3) {
		Vec3 p2 = m.nodes(c.node).pos;
		Vec3 p3 = m.nodes(m.corners(m.corners(c.next).opposite).node).pos;
		out = (5.0/12.0)*p1 - (1.0/12.0)*(p2+p3) + 0.75*p0;        
	} 
	else if(valence==4) {
		Vec3 p2 = m.nodes(m.corners(m.corners(c.next).opposite).node).pos;
		out = 0.375*p1 - 0.125*p2 + 0.75*p0;        
	}
	else {
		// rotate around extraordinary vertex,
		// calculate subdivision weights,
		// and interpolate vertex position
		double rv = 1.0/(double)valence;
		out = 0.0;
		int current = c.prev;
		for(int j=0; j<valence; j++) {
			double s = (0.25 + cos(2*M_PI*j*rv) + 0.5*cos(4*M_PI*j*rv))*rv;
			Vec3 p = m.nodes(m.corners(m.corners(current).prev).node).pos;

			out += s*p;
			current = m.corners(m.corners(m.corners(current).next).opposite).next;
		}
		out += 0.75* m.nodes(m.corners(c.prev).node).pos;        
	}
	return out;
}

// Modified Butterfly Subdivision Scheme from:
// Interpolating Subdivision for Meshes with Arbitrary Topology
// Denis Zorin, Peter Schroder, and Wim Sweldens
Vec3 ModifiedButterflySubdivision(Mesh& m, const Corner &ca, const Corner &cb, const Vec3& fallback)
{
	// calculate the valence of the two parent vertices
	int start = ca.prev;
	int current = start;
	int valenceA = 0;
	do {
		valenceA++;
		int op = m.corners(m.corners(current).next).opposite;
		if (op < 0) return fallback;
		current = m.corners(op).next;
	} 
	while(current != start);
	start = ca.next;
	current = start;
	int valenceB = 0;
	do {
		valenceB++;
		int op = m.corners(m.corners(current).next).opposite;
		if (op < 0) return fallback;
		current = m.corners(op).next;        
	} 
	while(current != start);

	// if both vertices have valence 6, use butterfly subdivision
	if(valenceA==6 && valenceB==6) {
		return ButterflySubdivision(m, ca,cb);
	}
	else if(valenceA==6)    // use a one-sided scheme
	{
		return OneSidedButterflySubdivision(m, valenceB,cb);
	}
	else if(valenceB==6)    // use a one-sided scheme
	{
		return OneSidedButterflySubdivision(m, valenceA,ca);
	}
	else    // average the results from two one-sided schemes
	{   
		return 0.5*(    OneSidedButterflySubdivision(m, valenceA,ca) 
					+   OneSidedButterflySubdivision(m, valenceB,cb)   );
	}
}

bool gAbort = false;

// collapse an edge on triangle "trinum".
// "which" is 0,1, or 2,
// where which==0 is the triangle edge from p0 to p1, 
// which==1 is the triangle edge from p1 to p2, 
// and which==2 is the triangle edge from p2 to p0, 
void CollapseEdge(Mesh& m, const int trinum, const int which, const Vec3 &edgevect, const Vec3 &endpoint,
				  vector<int> &deletedNodes, std::map<int,bool> &taintedTris, int &numCollapses, bool doTubeCutting)
{
	if (gAbort) return;
	// I wanted to draw a pretty picture of an edge collapse,
	// but I don't know how to make wacky angled lines in ASCII.
	// Instead, I will show the before case and tell you what needs to be done.

	//      BEFORE:
	//              *
	//             / \.
	//            /C0 \.
	//           /     \.
	//          /       \.
	//         /    B    \.
	//        /           \.
	//       /C1        C2 \.
	// P0 *---------------* P1
	//       \C2        C1 /
	//        \           /
	//         \    A    /
	//          \       /
	//           \     /
	//            \C0 /
	//             \ /
	//              *
	//
	//  We are going to collapse the edge between P0 and P1
	//      by deleting P1,
	//      and taking all references to P1,
	//          and rerouting them to P0 instead
	//
	//  What we need to do:
	//      Move position of P0
	//      Preserve connectivity in both triangles:
	//          (C1.opposite).opposite = C2.o
	//          (C2.opposite).opposite = C1.o
	//      Delete references to Corners of deleted triangles in both P0 and P1's Corner list
	//      Reassign references to P1:
	//          loop through P1 triangles: 
	//              rename P1 references to P0 in p lists.
	//              rename Corner.v references
	//      Copy P1's list of Corners over to P0's list of Corners
	//      Delete P1

	Corner ca_old[3], cb_old[3];
	ca_old[0] = m.corners(trinum, which);
	ca_old[1] = m.corners(ca_old[0].next);
	ca_old[2] = m.corners(ca_old[0].prev);
	bool haveB = false;
	if (ca_old[0].opposite>=0) {
		cb_old[0] = m.corners(ca_old[0].opposite);
		cb_old[1] = m.corners(cb_old[0].next);
		cb_old[2] = m.corners(cb_old[0].prev);
		haveB = true;
	}
	if (!haveB) {
		// for now, don't collapse
		return;
	}

	int P0 = ca_old[2].node;
	int P1 = ca_old[1].node;
	
	///////////////
	// avoid creating nonmanifold edges
	bool nonmanifold = false;
	bool nonmanifold2 = false;
	
	set<int>& ring0 = m.get1Ring(P0).nodes;
	set<int>& ring1 = m.get1Ring(P1).nodes;
	
	// check for intersections of the 1-rings of P0,P1
	int cl=0, commonVert=-1;
	for(set<int>::iterator it=ring1.begin(); it != ring1.end(); ++it)
		if (ring0.find(*it) != ring0.end()) {
			cl++;
			if (*it != ca_old[0].node && *it != cb_old[0].node) commonVert = *it;
		}
	
	nonmanifold = cl>2;
	nonmanifold2 = cl>3;
	
	if(nonmanifold && 
	   ca_old[1].opposite>=0 && cb_old[1].opposite>=0 &&
	   ca_old[2].opposite>=0 && cb_old[2].opposite>=0 ) // collapsing this edge would create a non-manifold edge
	{
		if(nonmanifold2) 
			return;
		
		bool topTet = false;
		bool botTet = false;
		// check if collapsing this edge will collapse a tet.
		if(m.corners(ca_old[1].opposite).node == m.corners(ca_old[2].opposite).node)
			botTet = true;
		
		if(m.corners(cb_old[1].opposite).node == m.corners(cb_old[2].opposite).node)   
			topTet = true;

		if(topTet^botTet) {
		
			// safe pyramid case.
			// collapse the whole tet!
			// First collapse the top of the pyramid,
			// then carry on collapsing the original verts.
			Corner cc_old[3],cd_old[3];
			if(botTet)  
				cc_old[0] = m.corners(ca_old[1].opposite);
			else // topTet
				cc_old[0] = cb_old[2];
			cc_old[1] = m.corners(cc_old[0].next);
			cc_old[2] = m.corners(cc_old[0].prev);
			if (cc_old[0].opposite<0) return;
			cd_old[0] = m.corners(cc_old[0].opposite);
			cd_old[1] = m.corners(cd_old[0].next);
			cd_old[2] = m.corners(cd_old[0].prev);
			int P2 = cc_old[2].node;
			int P3 = cc_old[1].node;
			
			// update tri props of all adjacent triangles of P0,P1 (do before CT updates!)
			for (int i=0; i<m.numTriChannels(); i++)                
				{};//TODO: handleTriPropertyEdgeCollapse(trinum, P2,P3,  cc_old[0], cd_old[0]);

			m.mergeNode(P2, P3);
			
			// Preserve connectivity in both triangles
			if (cc_old[1].opposite>=0)
				m.corners(cc_old[1].opposite).opposite = cc_old[2].opposite;
			if (cc_old[2].opposite>=0)
				m.corners(cc_old[2].opposite).opposite = cc_old[1].opposite;
			if (cd_old[1].opposite>=0)
				m.corners(cd_old[1].opposite).opposite = cd_old[2].opposite;
			if (cd_old[2].opposite>=0)
				m.corners(cd_old[2].opposite).opposite = cd_old[1].opposite;
			
			////////////////////
			// mark the two triangles and the one node for deletion
			int tmpTrinum = cc_old[0].tri;
			int tmpOthertri = cd_old[0].tri;
			m.removeTriFromLookup(tmpTrinum);
			m.removeTriFromLookup(tmpOthertri);
			taintedTris[tmpTrinum] = true;
			taintedTris[tmpOthertri] = true;
			deletedNodes.push_back(P3);

			numCollapses++;

			// recompute Corners for triangles A and B
			if(botTet)
				ca_old[0] = m.corners(ca_old[2].opposite);
			else
				ca_old[0] = m.corners(ca_old[1].prev);
			ca_old[1] = m.corners(ca_old[0].next);
			ca_old[2] = m.corners(ca_old[0].prev);
			cb_old[0] = m.corners(ca_old[0].opposite);
			cb_old[1] = m.corners(cb_old[0].next);
			cb_old[2] = m.corners(cb_old[0].prev);
		
			///////////////
			// avoid creating nonmanifold edges... again
			ring0 = m.get1Ring(ca_old[2].node).nodes;
			ring1 = m.get1Ring(ca_old[1].node).nodes;
			
			// check for intersections of the 1-rings of P0,P1
			cl=0;
			for(set<int>::iterator it=ring1.begin(); it != ring1.end(); ++it)
				if (*it != ca_old[0].node && ring0.find(*it) != ring0.end())
					cl++;
				
			if(cl>2) { // nonmanifold
				// this can happen if collapsing the first tet leads to another similar collapse that requires the collapse of a tet.
				// for now, just move on and pick this up later.

				// if the original component was very small, this first collapse could have led to a tiny piece of nonmanifold geometry.
				// in this case, just delete everything that remains.
				if(m.corners(ca_old[0].opposite).tri==cb_old[0].tri && m.corners(ca_old[1].opposite).tri==cb_old[0].tri && m.corners(ca_old[2].opposite).tri==cb_old[0].tri) {
					taintedTris[ca_old[0].tri] = true;
					taintedTris[cb_old[0].tri] = true;
					m.removeTriFromLookup(ca_old[0].tri);
					m.removeTriFromLookup(cb_old[0].tri);
					deletedNodes.push_back(ca_old[0].node);
					deletedNodes.push_back(ca_old[1].node);
					deletedNodes.push_back(ca_old[2].node);
				}
				return;
			}
		} else if(topTet && botTet && ca_old[1].opposite>=0 && ca_old[2].opposite>=0 && cb_old[1].opposite>=0 && cb_old[2].opposite>=0)
		{
			if(!(m.corners(ca_old[1].opposite).node == m.corners(ca_old[2].opposite).node && 
				m.corners(cb_old[1].opposite).node == m.corners(cb_old[2].opposite).node && 
				(m.corners(ca_old[1].opposite).node == m.corners(cb_old[1].opposite).node || 
				(m.corners(ca_old[1].opposite).node == cb_old[0].node && 
				m.corners(cb_old[1].opposite).node == ca_old[0].node) )))
			{
				// just collapse one for now.

				// collapse the whole tet!
				// First collapse the top of the pyramid,
				// then carry on collapsing the original verts.
				Corner cc_old[3],cd_old[3];
				
				// collapse top
				{
					cc_old[0] = m.corners(ca_old[1].opposite);
					cc_old[1] = m.corners(cc_old[0].next);
					cc_old[2] = m.corners(cc_old[0].prev);
					if (cc_old[0].opposite<0) return;
					cd_old[0] = m.corners(cc_old[0].opposite);
					cd_old[1] = m.corners(cd_old[0].next);
					cd_old[2] = m.corners(cd_old[0].prev);
					int P2 = cc_old[2].node;
					int P3 = cc_old[1].node;

					// update tri props of all adjacent triangles of P0,P1 (do before CT updates!)
					// TODO: handleTriPropertyEdgeCollapse(trinum, P2,P3,  cc_old[0], cd_old[0]);

					m.mergeNode(P2, P3);
					
					// Preserve connectivity in both triangles
					if (cc_old[1].opposite>=0)
						m.corners(cc_old[1].opposite).opposite = cc_old[2].opposite;
					if (cc_old[2].opposite>=0)
						m.corners(cc_old[2].opposite).opposite = cc_old[1].opposite;
					if (cd_old[1].opposite>=0)
						m.corners(cd_old[1].opposite).opposite = cd_old[2].opposite;
					if (cd_old[2].opposite>=0)
						m.corners(cd_old[2].opposite).opposite = cd_old[1].opposite;
					
					////////////////////
					// mark the two triangles and the one node for deletion
					int tmpTrinum = cc_old[0].tri;
					int tmpOthertri = cd_old[0].tri;
					taintedTris[tmpTrinum] = true;
					taintedTris[tmpOthertri] = true;
					m.removeTriFromLookup(tmpTrinum);
					m.removeTriFromLookup(tmpOthertri);
					deletedNodes.push_back(P3);

					numCollapses++;
				}
				// then collapse bottom
				{
					//cc_old[0] = [ca_old[1].opposite;
					cc_old[0] = cb_old[2];
					cc_old[1] = m.corners(cc_old[0].next);
					cc_old[2] = m.corners(cc_old[0].prev);
					if (cc_old[0].opposite<0) return;
					cd_old[0] = m.corners(cc_old[0].opposite);
					cd_old[1] = m.corners(cd_old[0].next);
					cd_old[2] = m.corners(cd_old[0].prev);
					int P2 = cc_old[2].node;
					int P3 = cc_old[1].node;

					// update tri props of all adjacent triangles of P0,P1 (do before CT updates!)
					// TODO: handleTriPropertyEdgeCollapse(trinum, P2,P3,  cc_old[0], cd_old[0]);
					
					m.mergeNode(P2, P3);

					// Preserve connectivity in both triangles
					if (cc_old[1].opposite>=0)
						m.corners(cc_old[1].opposite).opposite = cc_old[2].opposite;
					if (cc_old[2].opposite>=0)
						m.corners(cc_old[2].opposite).opposite = cc_old[1].opposite;
					if (cd_old[1].opposite>=0)
						m.corners(cd_old[1].opposite).opposite = cd_old[2].opposite;
					if (cd_old[2].opposite>=0)
						m.corners(cd_old[2].opposite).opposite = cd_old[1].opposite;
					
					////////////////////
					// mark the two triangles and the one node for deletion
					int tmpTrinum = cc_old[0].tri;
					int tmpOthertri = cd_old[0].tri;
					taintedTris[tmpTrinum] = true;
					taintedTris[tmpOthertri] = true;
					deletedNodes.push_back(P3);

					numCollapses++;
				}
				
				// Though we've collapsed a lot of stuff, we still haven't collapsed the original edge.
				// At this point we still haven't guaranteed that this original collapse weill be safe.
				// quit for now, and we'll catch the remaining short edges the next time this function is called.
				return;
			}
		}
		else if (doTubeCutting)
		{
			// tube case
			//cout<<"CollapseEdge:tube case" << endl;
			
			// find the edges that touch the common vert
			int P2 = commonVert;
			int P1P2=-1, P2P1, P2P0=-1, P0P2=-1; // corners across from the cutting seam
			int start = ca_old[0].next;
			int end = cb_old[0].prev;
			int current = start;
			do {
				// rotate around vertex P1 counter-clockwise
				int op = m.corners(m.corners(current).next).opposite;
				if (op < 0) errMsg("tube cutting failed, no opposite");
				current = m.corners(op).next;
				
				if(m.corners(m.corners(current).prev).node==commonVert)
					P1P2 = m.corners(current).next;
			}
			while(current != end);

			start = ca_old[0].prev;
			end = cb_old[0].next;
			current = start;
			do {
				// rotate around vertex P0 clockwise
				int op = m.corners(m.corners(current).prev).opposite;
				if (op < 0) errMsg("tube cutting failed, no opposite");
				
				current = m.corners(op).prev;
				if(m.corners(m.corners(current).next).node==commonVert)
					P2P0 = m.corners(current).prev;
			} while(current != end);

			if (P1P2 < 0 || P2P0 < 0) 
				errMsg("tube cutting failed, ill geometry");
			
			P2P1 = m.corners(P1P2).opposite;
			P0P2 = m.corners(P2P0).opposite;

			// duplicate vertices on the top half of the cut,
			// and use them to split the tube at this seam
			int P0b = m.addNode(Node(m.nodes(P0).pos));
			int P1b = m.addNode(Node(m.nodes(P1).pos));
			int P2b = m.addNode(Node(m.nodes(P2).pos));
			for (int i=0; i<m.numNodeChannels(); i++) {
				m.nodeChannel(i)->addInterpol(P0, P0, 0.5);
				m.nodeChannel(i)->addInterpol(P1, P1, 0.5);
				m.nodeChannel(i)->addInterpol(P2, P2, 0.5);
			}
			
			// offset the verts in the normal directions to avoid self intersections
			Vec3 offsetVec = cross(m.nodes(P1).pos-m.nodes(P0).pos, m.nodes(P2).pos-m.nodes(P0).pos);        
			normalize(offsetVec);
			offsetVec *= 0.01; // HACK: 
			m.nodes(P0).pos -= offsetVec;
			m.nodes(P1).pos -= offsetVec;
			m.nodes(P2).pos -= offsetVec;
			m.nodes(P0b).pos += offsetVec;
			m.nodes(P1b).pos += offsetVec;
			m.nodes(P2b).pos += offsetVec;

			// create a list of all triangles which touch P0, P1, and P2 from the top,
			map<int,bool> topTris;
			start = cb_old[0].next;
			end = m.corners(P0P2).prev;
			current = start;
			topTris[start/3]=true;
			do {
				// rotate around vertex P0 counter-clockwise
				current = m.corners(m.corners(m.corners(current).next).opposite).next;
				topTris[current/3]=true;
			} while(current != end);
			start = m.corners(P0P2).next;
			end = m.corners(P2P1).prev;
			current = start;
			topTris[start/3]=true;
			do {
				// rotate around vertex P0 counter-clockwise
				current = m.corners(m.corners(m.corners(current).next).opposite).next;
				topTris[current/3]=true;
			} while(current != end);
			start = m.corners(P2P1).next;
			end = cb_old[0].prev;
			current = start;
			topTris[start/3]=true;
			do {
				// rotate around vertex P0 counter-clockwise
				current = m.corners(m.corners(m.corners(current).next).opposite).next;
				topTris[current/3]=true;
			} while(current != end);
			
			// create two new triangles,
			int Ta = m.addTri(Triangle(P0,P1,P2));
			int Tb = m.addTri(Triangle(P1b,P0b,P2b));
			for (int i=0; i<m.numTriChannels(); i++) {
				m.triChannel(i)->addNew();
				m.triChannel(i)->addNew();
			}
					
			// sew the tris to close the cut on each side
			for(int c=0; c<3; c++) m.addCorner(Corner(Ta, m.tris(Ta).c[c]));
			for(int c=0; c<3; c++) m.addCorner(Corner(Tb, m.tris(Tb).c[c]));
			for(int c=0; c<3; c++) {
				m.corners(Ta,c).next = 3*Ta+((c+1)%3);
				m.corners(Ta,c).prev = 3*Ta+((c+2)%3);
				m.corners(Tb,c).next = 3*Tb+((c+1)%3);
				m.corners(Tb,c).prev = 3*Tb+((c+2)%3);
			}
			m.corners(Ta,0).opposite = P1P2;
			m.corners(Ta,1).opposite = P2P0;
			m.corners(Ta,2).opposite = ca_old[1].prev;
			m.corners(Tb,0).opposite = P0P2;
			m.corners(Tb,1).opposite = P2P1;
			m.corners(Tb,2).opposite = cb_old[1].prev;
			for (int c=0; c<3; c++) {
				m.corners(m.corners(Ta,c).opposite).opposite = 3*Ta+c;
				m.corners(m.corners(Tb,c).opposite).opposite = 3*Tb+c;
			}
			// replace P0,P1,P2 on the top with P0b,P1b,P2b.
			for(map<int,bool>::iterator tti=topTris.begin(); tti!=topTris.end(); tti++) {
				//cout << "H " << tti->first << " : " << m.tris(tti->first).c[0] << " " << m.tris(tti->first).c[1] << " " << m.tris(tti->first).c[2] << " " << endl;
				for(int i=0; i<3; i++) {
					int cn = m.tris(tti->first).c[i];
					set<int>& ring = m.get1Ring(cn).nodes;          
					
					if (ring.find(P0) != ring.end() && cn!=P0 && cn!=P1 && cn!=P2 && cn!=P0b && cn!=P1b && cn!=P2b) {
						ring.erase(P0);
						ring.insert(P0b);
						m.get1Ring(P0).nodes.erase(cn);
						m.get1Ring(P0b).nodes.insert(cn);
					}
					if (ring.find(P1) != ring.end() && cn!=P0 && cn!=P1 && cn!=P2 && cn!=P0b && cn!=P1b && cn!=P2b) {
						ring.erase(P1);
						ring.insert(P1b); 
						m.get1Ring(P1).nodes.erase(cn);
						m.get1Ring(P1b).nodes.insert(cn);                         
					}
					if (ring.find(P2) != ring.end() && cn!=P0 && cn!=P1 && cn!=P2 && cn!=P0b && cn!=P1b && cn!=P2b) {
						ring.erase(P2);
						ring.insert(P2b); 
						m.get1Ring(P2).nodes.erase(cn);
						m.get1Ring(P2b).nodes.insert(cn);                         
					}
					if(cn==P0) {
						m.tris(tti->first).c[i]=P0b;
						m.corners(tti->first,i).node = P0b;                        
						m.get1Ring(P0).tris.erase(tti->first);
						m.get1Ring(P0b).tris.insert(tti->first);
					}
					else if(cn==P1) {
						m.tris(tti->first).c[i]=P1b;
						m.corners(tti->first,i).node = P1b;
						m.get1Ring(P1).tris.erase(tti->first);
						m.get1Ring(P1b).tris.insert(tti->first);
					}
					else if(cn==P2) {
						m.tris(tti->first).c[i]=P2b;
						m.corners(tti->first,i).node = P2b;
						m.get1Ring(P2).tris.erase(tti->first);
						m.get1Ring(P2b).tris.insert(tti->first);
					}
				}
			}
			
			//m.sanityCheck(true, &deletedNodes, &taintedTris);
			
			return;
		}
		return;
	}
	if(ca_old[1].opposite>=0 && ca_old[2].opposite>=0 && cb_old[1].opposite>=0 && cb_old[2].opposite>=0 && ca_old[0].opposite>=0 && cb_old[0].opposite>=0 && 
		((m.corners(ca_old[1].opposite).node == m.corners(ca_old[2].opposite).node && // two-pyramid tubey case (6 tris, 5 verts)
		m.corners(cb_old[1].opposite).node == m.corners(cb_old[2].opposite).node && 
		(m.corners(ca_old[1].opposite).node == m.corners(cb_old[1].opposite).node || 
		(m.corners(ca_old[1].opposite).node==cb_old[0].node &&    // single tetrahedron case
		m.corners(cb_old[1].opposite).node==ca_old[0].node) )) 
		||
		(m.corners(ca_old[0].opposite).tri==m.corners(cb_old[0].opposite).tri && m.corners(ca_old[1].opposite).tri==m.corners(cb_old[0].opposite).tri && m.corners(ca_old[2].opposite).tri==m.corners(cb_old[0].opposite).tri   // nonmanifold: 2 tris, 3 verts
			&& m.corners(cb_old[0].opposite).tri==m.corners(ca_old[0].opposite).tri && m.corners(cb_old[1].opposite).tri==m.corners(ca_old[0].opposite).tri && m.corners(cb_old[2].opposite).tri==m.corners(ca_old[0].opposite).tri)
		))
	{
		// both top and bottom are closed pyramid caps, or it is a single tet
		// delete the whole component!
		// flood fill to mark all triangles in the component
		map<int,bool> markedTris;
		queue<int> triQ;
		triQ.push(trinum);
		markedTris[trinum] = true;
		int iters = 0;
		while(!triQ.empty()) {
			int trival = triQ.front();
			triQ.pop();
			for(int i=0; i<3; i++) {
				int newtri = m.corners(m.corners(trival,i).opposite).tri;
				if(markedTris.find(newtri)==markedTris.end()) {
					triQ.push(newtri);
					markedTris[newtri] = true;
				}
			}
			iters++;
		}
		map<int,bool> markedverts;
		for(map<int,bool>::iterator mit=markedTris.begin(); mit!=markedTris.end(); mit++) {
			taintedTris[mit->first] = true;
			markedverts[m.tris(mit->first).c[0]] = true;
			markedverts[m.tris(mit->first).c[1]] = true;
			markedverts[m.tris(mit->first).c[2]] = true;            
		}
		for(map<int,bool>::iterator mit=markedverts.begin(); mit!=markedverts.end(); mit++)
			deletedNodes.push_back(mit->first);
		return;
	}

	//////////////////////////
	// begin original edge collapse

	// update tri props of all adjacent triangles of P0,P1 (do before CT updates!)
	// TODO: handleTriPropertyEdgeCollapse(trinum, P0,P1,  ca_old[0], cb_old[0]);

	m.mergeNode(P0, P1);
	
	// Move position of P0
	m.nodes(P0).pos = endpoint + 0.5*edgevect;

	// Preserve connectivity in both triangles
	if (ca_old[1].opposite>=0)
		m.corners(ca_old[1].opposite).opposite = ca_old[2].opposite;
	if (ca_old[2].opposite>=0)
		m.corners(ca_old[2].opposite).opposite = ca_old[1].opposite;
	if (haveB && cb_old[1].opposite>=0)
		m.corners(cb_old[1].opposite).opposite = cb_old[2].opposite;
	if (haveB && cb_old[2].opposite>=0)
		m.corners(cb_old[2].opposite).opposite = cb_old[1].opposite;
	
	////////////////////
	// mark the two triangles and the one node for deletion
	taintedTris[ca_old[0].tri] = true;
	m.removeTriFromLookup(ca_old[0].tri);
	if (haveB) {
		taintedTris[cb_old[0].tri] = true;
		m.removeTriFromLookup(cb_old[0].tri);    
	}
	deletedNodes.push_back(P1);
	numCollapses++;
}

} // namespace
